home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / DCLAP 6d / dclap6d / SeqPups / appsrc / autoseq.src / CPeakList.C < prev    next >
Text File  |  1996-07-05  |  21KB  |  725 lines
  1. //    ============================================================================
  2. //    CPeakList.C                                                        80 columns
  3. //    Reece Hart    (reece@ibc.wustl.edu)                                tab=4 spaces
  4. //    Washington University School of Medicine, St. Louis, Missouri
  5. //    This source is hereby released to the public domain.  Bug reports, code
  6. //    contributions, and suggestions are appreciated (to email address above).
  7. //    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -
  8. //    Please see CPeakList.H for a description of the CPeakList class.
  9. //    ========================================|===================================
  10.  
  11.  
  12. #include    "CPeakList.H"
  13. #ifndef WIN_MAC
  14. #include    "CTrace.H"
  15. #else
  16. #include    "CTrace.C"  // dgg
  17. #endif
  18. #include    "RInlines.H"
  19. #include    <math.h>
  20. #include    <stdlib.h>
  21. #include    <stdio.h>
  22.  
  23. vrsn CPeakList::version = "94.05.03";
  24.  
  25. CPeakList::CPeakList(void):
  26.     hmean(0),
  27.     hvariance(0),
  28.  
  29.     hm0(0),
  30.     hmdecay(0),
  31.  
  32.     hv0(0),
  33.     hvdecay(0),
  34.  
  35.     w0(0),
  36.     wconst(0),
  37.  
  38.     group1_left(0),
  39.     group1_leftidx(0),
  40.     group1_right(0),
  41.     group1_rightidx(0),
  42.     group1_hmean(0),
  43.     group1_hvariance(0),
  44.     group1_wd_mean(0),
  45.     group1_index_mean(0),
  46.  
  47.     group2_left(0),
  48.     group2_leftidx(0),
  49.     group2_right(0),
  50.     group2_rightidx(0),
  51.     group2_hmean(0),
  52.     group2_wd_mean(0),
  53.     group2_hvariance(0),
  54.     group2_index_mean(0),
  55.  
  56.     FTrace(NULL),
  57.     PTrace(NULL),
  58.     RTrace(NULL)
  59.     {}
  60.  
  61. CPeakList::~CPeakList(void)
  62.     {
  63.     // no error to delete NULL
  64.     delete FTrace;
  65.     delete PTrace;
  66.     delete RTrace;
  67.     }
  68.  
  69.  
  70. //    ============================================================================
  71. //    Analyze
  72. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  73. //    Automates the process of computing peak probabilities from a peak list.
  74. //    ========================================|===================================
  75. bool
  76. CPeakList::Analyze(
  77.     CTrace<trace_t>& ct)
  78.     {
  79.     if (ct.Peaks().Size() != 0)
  80.         {
  81.         ComputeFTrace(ct.Size());
  82.         FTrace->LeftCutoff(ct.LeftCutoff());
  83.         FTrace->RightCutoff(ct.RightCutoff());
  84.  
  85.         ComputeRTrace(ct);
  86.         RTrace->LeftCutoff(ct.LeftCutoff());
  87.         RTrace->RightCutoff(ct.RightCutoff());
  88.  
  89.         CalculateStats();
  90.         CalculatePeakStats(ct.Min());
  91.         }
  92.     return TRUE;
  93.     }
  94.  
  95. //    ============================================================================
  96. //    CalculateStats
  97. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  98. //    This is a big ugly function.  What we need to determine is:
  99. //        1. the overall mean and variance of peak heights
  100. //        2. the linear increase of peak width with index
  101. //        3. the exponential decay equation for the mean (& variance) of peak
  102. //             heights
  103. //    For 1, we just iterate over the whole list and get the mean & var.  Easy.
  104. //    For 2 & 3, we wanted to divide the peak list population into 3 groups, and
  105. //    fit the equations to means from the first and last third.  This turns out
  106. //    to be very inaccurate because of large fluctuations in both height and
  107. //    width in the last third.  We therefore resort to fitting the equation
  108. //    between the first and second thirds of the population.
  109. //    We initially tried to use the peak widths to model the decay with an
  110. //    assumption of constant peak area, but for the same reason, this method was
  111. //    error prone.
  112. //    ========================================|===================================
  113. void
  114. CPeakList::CalculateStats(void)
  115.     {
  116.     uint        index;
  117.     uint        population_size;
  118.     double        sum;
  119.     CPeakList&    peaks = *this;
  120.  
  121.     if (size < 2)
  122.         return;
  123.  
  124.     //
  125.     // compute the mean & variance of height for all peaks (not groups)
  126.     //
  127.     for (sum=0,index=0; index < size; index++)
  128.         sum += peaks[index].height;
  129.     hmean = sum/size;
  130.     for (sum=0,index=0; index < size; index++)
  131.         sum += pow(peaks[index].height - hmean,2);
  132.     hvariance = sum/(size-1);
  133.  
  134.  
  135.     //
  136.     // divide peaks into two groups
  137.     // the groups will be considered as two populations for the purpose
  138.     // of calculating means, etc.  From those, we'll try to fit an
  139.     // exponential decay to the height and a linear function to width
  140.     // (by least squares).
  141.     //
  142.     group1_left        = 0;                    // group 1 = first third
  143.     group1_right    = (uint)size/3-1;
  144.     group1_leftidx    = (*this)[group1_left].center;
  145.     group1_rightidx    = (*this)[group1_right].center;
  146.  
  147.     group2_left        = group1_right + 1;        // group 2 = second third
  148.     group2_right    = (uint)size*2/3;
  149.     group2_leftidx    = (*this)[group2_left].center;
  150.     group2_rightidx    = (*this)[group2_right].center;
  151.  
  152.  
  153.     //
  154.     // model peak width expansion by linear least squares fit of widths
  155.     // ie.  w(i) = wd0 + wd_const * i;
  156.     //
  157.     double    si    =0;                            // sum of index
  158.     double    si2    =0;                            // sum of index^2
  159.     double    sw    =0;                            // sum of width
  160.     double    sw2    =0;                            // sum of width^2
  161.     double    siw    =0;                            // sum of width * index
  162.     for (index=0; index < size; index++)
  163.         {
  164.         si    += peaks[index].center;
  165.         si2    += pow(peaks[index].center,2);
  166.         sw    += peaks[index].width;
  167.         siw    += peaks[index].width * peaks[index].center;
  168.         }
  169.     wconst = (si*sw - size*siw)/(pow(si,2)*si2);
  170.     w0 = (siw - si2*wconst)/si;
  171.  
  172.  
  173.     //
  174.     // model peak height decay with exponential
  175.     // ie.  corrected height = hc(i) = hm0 * exp(-h_decay*i)
  176.     // 1. divide peaks into two distinct groups (as above)
  177.     // 2. compute mean height & variance and mean index for each group
  178.     //        (index is necessary to compensate for skewed distribution of peaks)
  179.     // 3. fit exponential decay to mean from the two groups
  180.     // 4. compute hm0, the height at i=0
  181.     //
  182.  
  183.     //
  184.     // group 1 calculations
  185.     //
  186.     population_size = group1_right - group1_left + 1;
  187.  
  188.     // height mean
  189.     for (sum=0,index=group1_left; index <= group1_right; index++)
  190.         sum += peaks[index].height;
  191.     group1_hmean = sum/population_size;
  192.  
  193.     // height variance
  194.     for (sum=0,index=group1_left; index <= group1_right; index++)
  195.         sum += pow(peaks[index].height - group1_hmean,2);
  196.     group1_hvariance = sum/(population_size-1);
  197.  
  198.     // width mean
  199.     for (sum=0,index=group1_left; index <= group1_right; index++)
  200.         sum += peaks[index].width;
  201.     group1_wd_mean = sum/population_size;
  202.  
  203.     // index mean
  204.     for (sum=0,index=group1_left; index <= group1_right; index++)
  205.         sum += peaks[index].center;
  206.     group1_index_mean = sum/population_size;
  207.  
  208.  
  209.     //
  210.     // group 2 calculations
  211.     //
  212.     population_size = group2_right - group2_left + 1;
  213.  
  214.     // height mean
  215.     for (sum=0,index=group2_left; index <= group2_right; index++)
  216.         sum += peaks[index].height;
  217.     group2_hmean = sum/population_size;
  218.  
  219.     // height variance
  220.     for (sum=0,index=group2_left; index <= group2_right; index++)
  221.         sum += pow(peaks[index].height - group2_hmean,2);
  222.     group2_hvariance = sum/(population_size-1);
  223.  
  224.     // width mean
  225.     for (sum=0,index=group2_left; index <= group2_right; index++)
  226.         sum += peaks[index].width;
  227.     group2_wd_mean = sum/population_size;
  228.  
  229.     // index mean
  230.     for (sum=0,index=group2_left; index <= group2_right; index++)
  231.         sum += peaks[index].center;
  232.     group2_index_mean = sum/population_size;
  233.  
  234.  
  235.     //
  236.     // fit exponential
  237.     // determine hmdecay, hm0, hvdecay, and hv0
  238.     //
  239.     // the following attempts to use peak width to preserve area under peak
  240.     // hdecay = log(group1_hmean*group1_wd_mean / group2_hmean*group2_wd_mean)/
  241.     //            (group2_index_mean - group1_index_mean);
  242.  
  243.     hmdecay = log(group1_hmean/group2_hmean)/
  244.                             (group2_index_mean - group1_index_mean);
  245.     hm0 = group1_hmean/exp(-hmdecay * group1_index_mean);
  246.     hvdecay = log(group1_hvariance/group2_hvariance)/
  247.                             (group2_index_mean - group1_index_mean);
  248.     hv0 = group1_hvariance/exp(-hvdecay * group1_index_mean);
  249.     }
  250.  
  251.  
  252. //    ============================================================================
  253. //    CalculatePeakStats
  254. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  255. //    Computes P(H), P(D|NULL), P(D|H), and P(H|D) for all peaks.
  256. //
  257. //    Peaks are initially chosen by a concave down method (see PickPeaks).  From
  258. //    these initial selections, the mean and variance of peak heights are
  259. //    calculated, as are the terms for an exponential decay of peak height with
  260. //    trace position and terms for a linear expansion of peak width.
  261. //
  262. //    h            = arbitrary height (ie. h=h(3)=height at position 3)
  263. //    h(i)        = height at index i
  264. //    h_mean        = mean of peak heights computed by CalculateStats
  265. //    h_mean(i)    = mean of peak heights corrected for exponential decay
  266. //    h_var        = variance of peak heights computed by CalculateStats
  267. //    h_baseline    = baseline height
  268. //    noise_var    = variance in signal noise; determined empirically by
  269. //                    subtracting the FTrace from observed signal.
  270. //    FT(i)        = expected fluorescence trace of chosen peaks (see
  271. //                    ComputeFTrace)
  272. //
  273. //    P(H), P(D|H), and P(D|NULL) are computed using some method (see below) and
  274. //    then plugged into Bayes' Theorem.
  275. //      P(H|D) = P(H) * P(D|H) / P(D|NULL)
  276. //    ========================================|===================================
  277. void
  278. CPeakList::CalculatePeakStats(
  279.     double    h_baseline)
  280.     {
  281.     // we need the FTrace and RTrace; exit if they're not NULL
  282.     if ((NULL == FTrace) || (NULL == RTrace))
  283.         return;
  284.  
  285.     double noisevar=RTrace->Variance();
  286.  
  287.     for (uint index=0;index<size;index++)
  288.         {
  289.         PeakRec& peak = (*this)[index];        // current peak record
  290.  
  291.         // height & height variance decay exponentially
  292.         double    fluorexp    = ExpDecay(hm0,hmdecay,index);
  293.         double    fluorvar    = ExpDecay(hv0,hvdecay,index);
  294.  
  295.         double    varprod  = hvariance * noisevar;
  296.  
  297.         // the following are terms defined in the project report and have
  298.         // no particular meaning other than to aggregate terms.
  299.         // See the project report (p. 27) for details.
  300.         double    A,B,C,D,E,F;
  301.         A = (hvariance + noisevar)/2/varprod;
  302.         B = (hvariance*(peak.height-h_baseline) + noisevar*fluorexp)/varprod;
  303.         C = -(hvariance*pow((peak.height-h_baseline),2) + 
  304.                         noisevar*pow(fluorexp,2)) / 2 / varprod;
  305.         D = sqrt(A);
  306.         E = B/2/D;
  307.         F = C-pow(B,2)/4/A;
  308.  
  309.         peak.PH  = 1.0;
  310. //        peak.PH  = GaussianV((peak.height-h_baseline),fluorexp,fluorvar);
  311.         peak.PDN = GaussianV(peak.height-h_baseline,0,noisevar);
  312. //        peak.PDH = GaussianV(peak.height-h_baseline-fluorexp,0,noisevar);
  313.  
  314. #ifdef WIN_MAC
  315.         peak.PDH = exp(-F)/4/D/sqrt(varprod*_PI);  // dgg - no erfc() ..
  316. #else
  317.         peak.PDH = exp(-F)/4/D/sqrt(varprod*_PI)*erfc(E);   
  318. #endif
  319. // debug print: printf("%f\t%f\t%f\t%f\t%f\t%f\n",A,B,C,D,E,F);
  320.  
  321.         // and finally Bayes' Theorem:
  322. //        peak.PHD = peak.PH * peak.PDH / peak.PDN;
  323.         peak.PHD = peak.PH * peak.PDH / ( peak.PDH+peak.PDN);
  324.         }
  325.     }
  326.  
  327. //    ============================================================================
  328. //    WriteDeltas
  329. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  330. //    For every subsequence of length nbhd, write the index of the 3' most peak
  331. //    center for that subsequence, the subsequence itself, and the delta for the
  332. //    subsequence.
  333. //    The delta between bases m and n is defined to be the distance between the
  334. //    centers of their corresponding peaks m and n in the trace.
  335. //    ========================================|===================================
  336. void
  337. CPeakList::WriteDeltas(
  338.     ostream& os,
  339.     uint    nbhd,
  340.     bool    header)
  341.     {
  342.     char*    seq = Sequence();
  343.  
  344.     if (header)
  345.         os    << "Pk Ctr"
  346.             << "\t"
  347.             << "Delta"
  348.             << "\t"
  349.             << "Sequence"
  350.             << endl;
  351.  
  352.     for(ulong index=nbhd;index<size;index++)
  353.         {
  354.         os    << (*this)[index].center
  355.             << "\t"
  356.             << Delta(index,nbhd)
  357.             << "\t";
  358.         for (ulong base=index-nbhd+1;base<=index;base++)
  359.             os << seq[base];
  360.         os    << endl;
  361.         }
  362.     delete seq;
  363.     }
  364.  
  365.  
  366. //    ============================================================================
  367. //    ComputeFTrace
  368. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  369. //    Computes a trace which represents the expected fluorescence from the list
  370. //    of peaks by assuming a Gaussian distribution at each peak (using the center,
  371. //    height, and width fields of the PeakRec).  The extent parameter specifies
  372. //    the horizontal extent of the Gaussian on each side of center.  The size of
  373. //    the original trace must be passed so that the original and fluorescence
  374. //    traces will be the same size.
  375. //    ========================================|===================================
  376. inline
  377. double
  378. PeakModel(
  379.     trace_t        height,
  380.     tracepos    center,
  381.     tracepos    width,
  382.     tracepos    index)
  383.     {
  384.     //                                  2
  385.     //                      -(i-center)  / (2 x width)
  386.     // F(i,p) =    height x e
  387.     //
  388.     // F(i,p) is fluorescence at index i as a result of peak p
  389.     // height = height of peak p
  390.     // center = center of peak p
  391.     // width  = width of peak p
  392.     // 
  393.     return height * exp(-pow((index-center),2)/width/2);
  394.     }
  395.  
  396. CTrace<double>*
  397. CPeakList::ComputeFTrace(
  398.     size_t        sz,
  399.     tracepos    extent)
  400.     {
  401.     delete FTrace; FTrace = NULL;            // recreate from scratch
  402.  
  403.     if (FTrace = new CTrace<double>(sz))
  404.         {
  405.         CTrace<double>&    trace=*FTrace;
  406.         tracepos index;
  407.         for (index=0; index<sz; index++)    // zero all indices of the FT
  408.             trace[index]=0;
  409.  
  410.         for (size_t peak_no=0; peak_no<size; peak_no++)
  411.             {
  412.             PeakRec&    peak = (*this)[peak_no];    // the current peak to
  413.                                                     // model with Gaussian
  414.  
  415.             tracepos leftBound = peak.center-extent;
  416.             tracepos rightBound = peak.center+extent;
  417.  
  418.             if (leftBound<0)        leftBound = 0;
  419.             if (rightBound>sz-1)    rightBound = sz-1;
  420.  
  421.             for (index=leftBound;index<=rightBound;index++)
  422.                 trace[index] +=
  423.                         PeakModel(peak.height,peak.center,peak.width,index);
  424.             }
  425.  
  426.         FTrace->CalculateStats();
  427.         }
  428.     return FTrace;
  429.     }
  430.  
  431.  
  432. //    ============================================================================
  433. //    ComputePTrace
  434. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  435. //    Generates a qualitative trace which shows peak height and width for every
  436. //    peak in the peak list. The argument passed to ComputePTrace is the size of
  437. //    the original trace, which will also be the size of the peak trace. 
  438. //    ========================================|===================================
  439. CTrace<double>*
  440. CPeakList::ComputePTrace(
  441.     CTrace<trace_t>&    src_trace)
  442. //    size_t        sz)
  443.     {
  444.     delete PTrace; PTrace = NULL;            // recreate from scratch
  445.  
  446.     tracepos    sz=src_trace.Size();
  447.     tracepos    offset=src_trace.LeftCutoff();
  448.  
  449.     if (PTrace = new CTrace<double>((size_t)sz))
  450.         {
  451.         CTrace<double>& ptrace = *PTrace;
  452.  
  453.         for (tracepos index = 0; index<sz; index++)
  454.             ptrace[index]=0;
  455.  
  456.         for (uint peak_num=0; peak_num < size; peak_num++)
  457.             {
  458.             PeakRec&    pr = (*this)[peak_num];
  459.             double        height = pr.height;
  460.             double        half_height = height/2.0;
  461.             tracepos    midpoint = pr.center;
  462.             tracepos    left = (tracepos)pr.leftBound;
  463.             tracepos    right = (tracepos)pr.rightBound;
  464.         
  465.             for (tracepos index=left; index<=right; index++)
  466.                 if (index<sz)    ptrace[index+offset]=half_height;
  467.  
  468.             if (midpoint<sz)    ptrace[midpoint+offset]=height;
  469.             }
  470.  
  471.         PTrace->LeftCutoff(src_trace.LeftCutoff());
  472.         PTrace->RightCutoff(src_trace.RightCutoff());
  473.  
  474.         PTrace->CalculateStats();
  475.         }
  476.  
  477.     return PTrace;
  478.     }
  479.  
  480. //    ============================================================================
  481. //    ComputeRTrace
  482. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  483. //    Generates a new trace of residuals (data that cannot be explained by peaks)
  484. //    by subtracting the computed FTrace from the original trace.  A pointer to
  485. //    the result is stored in RTrace and returned (or NULL if error).
  486. //    ========================================|===================================
  487. CTrace<double>*
  488. CPeakList::ComputeRTrace(
  489.     CTrace<trace_t>& src_trace)
  490.     {
  491.     // FTrace must already be computed
  492.     if (NULL == FTrace)
  493.         return NULL;
  494.  
  495.     delete PTrace; PTrace = NULL;            // recreate from scratch
  496.  
  497.     if (RTrace = new CTrace<double>(src_trace.Size()))
  498.         {
  499.         for (uint index=0;index<src_trace.Size();index++)
  500.             (*RTrace)[index] = src_trace[index]-(*FTrace)[index];
  501.     
  502.         RTrace->CalculateStats();
  503.         }
  504.  
  505.     return RTrace;        
  506.     }
  507.  
  508.  
  509. //    ============================================================================
  510. //    Offset
  511. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  512. //    Translate the trace position-dependent fields of every PeakRec to account
  513. //    for the left cutoff.
  514. //    ========================================|===================================
  515. void
  516. CPeakList::Offset(
  517.     tracepos    offset)
  518.     {
  519.     for (uint index = 0; index < size ; index++)
  520.         {
  521.         PeakRec&    pr = (*this)[index];
  522.         pr.center += offset;
  523.         pr.leftBound += offset;
  524.         pr.rightBound += offset;
  525.         }
  526.  
  527.     group1_leftidx += offset;
  528.     group1_rightidx += offset;
  529.     group1_index_mean += offset;
  530.     group2_leftidx += offset;
  531.     group2_rightidx += offset;
  532.     group2_index_mean += offset;
  533.     }    
  534.  
  535.  
  536. //    ============================================================================
  537. //    Sequence
  538. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  539. //    Writes the NULL-terminated sequence to buf.  If buf is NULL on entry, it
  540. //    will be created and it is the caller's responsibility to deallocate it.
  541. //    If successful, a pointer to buf is returned; otherwise NULL is returned.
  542. //    NOTE: use delete to deallocate, NOT free.
  543. //    ========================================|===================================
  544. char*
  545. CPeakList::Sequence(char* buf)
  546.     {
  547.     if (buf == NULL)
  548.         buf=new char[size+1];
  549.  
  550.     if (buf != NULL)
  551.         {
  552.         uint i;
  553.         for (i = 0; i<size; i++)
  554.             buf[i]=DNA_BASES[(*this)[i].whichTrace];
  555.         buf[i]=NULL;                        // i = size (not size-1) upon exit
  556.         }
  557.  
  558.     return buf;
  559.     }
  560.  
  561. //    ============================================================================
  562. //    operator<<
  563. //    -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -
  564. //    Display a summary of the peak list and it's members
  565. //    ========================================|===================================
  566. ostream&
  567. operator<<(
  568.     ostream&    os,
  569.     CPeakList&    cpl)
  570.     {
  571.     os    << "Peak report" << endl;
  572.     os    << "-----------" << endl;
  573.  
  574.     if (cpl.size == 0)
  575.         return os << "None picked." << endl;
  576.  
  577.     // else...
  578.     os    << "Number picked:\t" << cpl.size << endl
  579.         << "height:\t\t" << cpl.hmean << "+/-" << sqrt(cpl.hvariance) << endl
  580.         << "ht mean decay:\t" << cpl.hm0 << " * exp(-"
  581.                                 << cpl.hmdecay << " * x)" << endl
  582.         << "ht var decay:\t" << cpl.hv0 << " * exp(-"
  583.                                 << cpl.hvdecay << " * x)" << endl
  584.         << "width:\t\t" << cpl.w0 << " + " << cpl.wconst << " * x" << endl;
  585.     os    << "group1:" << endl
  586.         << "  peaks:      [" << cpl.group1_left << "," << cpl.group1_right
  587.             << "]" << endl
  588.         << "  indices:   [" << cpl.group1_leftidx << ","
  589.                              << cpl.group1_rightidx << "]" << endl
  590.         << "  height =     " << cpl.group1_hmean << "+/-"
  591.                              << sqrt(cpl.group1_hvariance) << endl
  592.         << "  wdmean =     " << cpl.group1_wd_mean << endl
  593.         << "  index mean = " << cpl.group1_index_mean << endl;
  594.     os    << "group2:" << endl
  595.         << "  peaks:      [" << cpl.group2_left << "," << cpl.group2_right
  596.             << "]" << endl
  597.         << "  indices:   [" << cpl.group2_leftidx << ","
  598.                              << cpl.group2_rightidx << "]" << endl
  599.         << "  height =     " << cpl.group2_hmean << "+/-"
  600.                              << sqrt(cpl.group2_hvariance) << endl
  601.         << "  wdmean =     " << cpl.group2_wd_mean << endl
  602.         << "  index mean = " << cpl.group2_index_mean << endl;
  603.  
  604.     if (cpl.RTrace)
  605.         os << "noise:\t\t" << cpl.RTrace->Mean()
  606.             << "+/-" << sqrt(cpl.RTrace->Variance()) << endl;
  607.  
  608.     if (cpl.size != 0)
  609.         {
  610.         os << setw(3) << "#" << "  " << PeakRecHeader;
  611.         uint i=1;
  612.         for (SeqNode<PeakRec>* sn = cpl.first;
  613.             NULL != sn;
  614.             sn=sn->next, i++)
  615.             os << setw(3) << i << "  " << sn->data;
  616.  
  617.         return os;
  618.         }
  619.  
  620.     return os;
  621.     }
  622.  
  623.  
  624. //    ============================================================================
  625. //    PeakRec::PeakRec (constructor)
  626. //    ========================================|===================================
  627. PeakRec::PeakRec(void):
  628.         whichTrace(X),
  629.         center(0),
  630.         height(0),
  631.         width(0),
  632.         leftBound(0),
  633.         rightBound(0),
  634.         PH(0),
  635.         PDN(0),
  636.         PDH(0),
  637.         PHD(0)
  638.         {}
  639.  
  640.  
  641. //    ============================================================================
  642. //    PeakRecHeader Manipulator and operator<< for PeakRec
  643. //    ========================================|===================================
  644. const int BASE_FW=6;                        // base field width
  645. const int INDEX_FW=8;                        // index field width
  646. const int PROB_FW=10;                        // width for P(X) fields
  647.  
  648. ostream&
  649. PeakRecHeader(
  650.     ostream& os)
  651.     // ostream manipulator to put column headings onto stream
  652.     {
  653.     char    buf[120];
  654.     sprintf(buf,"%-2s%5s%7s%8s%8s%7s%10s%10s%10s%10s\n",
  655.         "X",
  656.         "ht",
  657.         "ctr",
  658.         "left",
  659.         "right",
  660.         "width",
  661.         "P(H)",
  662.         "P(D|NULL)",
  663.         "P(D|H)",
  664.         "P(H|D)");
  665.     return os << buf;
  666. #if (FALSE)
  667.     return os
  668.         << setiosflags(ios::left)
  669.         << setw(BASE_FW)    << "base"
  670.         << setiosflags(ios::right)
  671.         << setw(INDEX_FW)    << "height"
  672.         << setw(INDEX_FW)    << "center"
  673.         << setw(INDEX_FW)    << "left"
  674.         << setw(INDEX_FW)    << "right"
  675.         << setw(INDEX_FW)    << "width"
  676.         << setw(PROB_FW)    << "P(H)"
  677.         << setw(PROB_FW)    << "P(D|NULL)"
  678.         << setw(PROB_FW)    << "P(D|H)"
  679.         << setw(PROB_FW)    << "P(H|D)"
  680.         << endl;
  681. #endif
  682.     }
  683.  
  684. ostream&
  685. operator<<(
  686.     ostream& os,
  687.     const PeakRec& pr)
  688.     {
  689.     char    buf[120];
  690.     sprintf(buf,"%-2c%5.0f%7d%8.1f%8.1f%7.1f%10.2e%10.2e%10.2e%10.2e\n",
  691.         DNA_BASES[pr.whichTrace],
  692.         pr.height,
  693.         pr.center,
  694.         pr.leftBound,
  695.         pr.rightBound,
  696.         pr.width,
  697.         pr.PH,
  698.         pr.PDN,
  699.         pr.PDH,
  700.         pr.PHD);
  701.     return os << buf;
  702. #if (FALSE)
  703.     return os
  704.         << setiosflags(ios::right)
  705.         << setw(BASE_FW)    << DNA_BASES[pr.whichTrace]
  706.         << setprecision(5)
  707.         << resetiosflags(ios::fixed|ios::scientific)
  708.         << setiosflags(ios::right)
  709.         << setw(INDEX_FW)    << pr.height
  710.         << setw(INDEX_FW)    << pr.center
  711.         << setw(INDEX_FW)    << pr.leftBound
  712.         << setw(INDEX_FW)    << pr.rightBound
  713.         << setw(INDEX_FW)    << pr.width
  714.         << setprecision(2)                    // always use 2-place scientific
  715.         << setiosflags(ios::scientific)        // notation for probabilities
  716.         << setw(PROB_FW)    << pr.PH
  717.         << setw(PROB_FW)    << pr.PDN
  718.         << setw(PROB_FW)    << pr.PDH
  719.         << setw(PROB_FW)    << pr.PHD
  720.         << endl;
  721. #endif
  722.     }
  723.  
  724.  
  725.